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Abstract.  The  paper  addresses  the  problem  of  plate  modeling  within  the 
framework  of  the  h-p  version  of  the  finite  element  method.  A  natural 
hierarchy  of  models  is  constructed.  The  lowest  member  of  the  hierarchy  is 
the  well  known  Reissner-Mindl in  model.  It  is  shown  that  higher  degree 
elements  do  not  show  any  locking  effects  for  the  models  under  consideration. 

Key  words:  Plates,  Reissner-Mindl in  model,  h-p  version  of  finite  element 
method 

1 •  Introduction. 

The  h-p  version  of  the  finite  element  method  was  developed  during  the 
last  10  years.  The  h-p  version  programs  MSC/PROBE,  FIESTA  are  on  the  market 
and  the  h-p  version  research  code  STRIPE  (Aeronautical  Institute  of  Sweden) 
is  used  in  industry  as  well. 

The  h-p  code  provides  large  flexibility.  It  allows  hierarchical 
modeling  of  the  plates  and  shells,  avoids  the  locking  problems  and  provides 
tools  for  effective  quality  control. 

This  paper  addresses  a  natural  hierarchical  modeling  which  the  h-p 
version  provides,  including  the  Reissner-Mindl in  model.  In  addition  it  shows 
that  the  h-p  version  avoids  the  problem  of  locking  which  arises  in  the 
standard  finite  element  method. 


2.  The  plate  problem  and  its  hierarchical  mode  1 s . 

3  d  d 

Let  fl  =  (x  =  (Xj,  x2>  x3)  €  R  ,  (x1,x2)  c  u f  ~2  <  x3  <  2  }  be  the 

three  dimensional  plate  w  with  the  thickness  d.  Further  we  let  S  =  {x  € 

3  d  d  3  d 

IR  ,  (Xj,x2)  e  T,  -2  <  x3  <  2  >  and  R+  *  {  x  e  R  ,  (x^.Xg)  e  u,  x3  =  ±2  } 
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where  r  denotes  the  boundary  of  w. 
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We  will  consider  the  three  dimensional  elasticity  problem  on  fl  with 
the  Hooke’s  Law 
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where  <rij,Tij  sure  the  normal  and  tangential  stresses, 
are  the  strain  components,  and  u^,  1  =  1,2,3  u  = 
the  displacements. 

For  isotropic  materials  we  have  A  =  {a^>,  1 , J  = 


respectively,  e^, 
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=  M  = 


2(1  +  v)  ’ 


A.p  are  the  Lam6  constants,  E  is  the  modulus  of  elasticity  and  v  is  the 
Poisson  ratio.  As  usual,  the  strain  energy  is 


(2.2) 


J 

n 


where  the  stress-strain  relations  in  (2.2)  are  governed  by  (2.1)  with  the 
matrix  A  =  {a^}.  A  is  used  as  an  upper  index  in  the  expression  (2.2)  to 
emphasize  this  dependence. 

We  will  assume  that  the  plate  is  homogeneous,  i.e.  that  are 

constants.  Further  we  will  assume  that  half  of  the  load  qtx^.x^)  is  on  R+ 
and  half  on  R  .  The  total  energy  is  then 
(2.3)  GA(u)  =  eA(u)  -  Q(u) 


2 


where 


(2.5)  Q(u)  =  Jl[u3(x1*x2*i)  *  u3(xi‘x2*"l)]dx 

u 

The  exact  solution  u  =  of  the  plate  problem  Is  the 

A  13 

mlnimizer  of  G  (u)  over  the  subspace  H(Q)  c  (H  (Q) )  where  K(H) 

1  3 

constraints  (H  (fl))  on  S  (not  on  R+).  The  boundary  conditions  of  the 
plate  problem  are  uniquely  characterized  by  H(G)  (respectively  the 
constraints  of  (H*(n))3  on  S). 

Let  now  n  =  (n^n^Oj),  ^  >  0  Integer.  Then  by  the  solution  of  the 

a"’) 

g 

G  (u)  over  the  subspace  Jf(n)  c  i'(£l)  of  all  functions  of  the  form 

ni  J 

(2.6)  u[n)(x)  =  £  ui,j(xrx2)[¥]  ’  1  =  1*2’3 

J=0 

and  the  matrix  B  =  {b^>  Is  used  as  Hooke’s  law.  In  general,  B  *  A. 
Because  of  our  symmetry  assumptions  we  made  earlier  we  have 

Uj  =  u2  j  =  0  for  J  even 

and  u_  =  0  for  j  odd. 

J*  J 

B  ( n ) 

The  solution  u  ,  Is  the  approximation  of  the  exact  solution  u  of  the 
three  dimensional  problem. 

We  can  use  various  n  for  the  approximate  solution  and  study  the  rate 
of  convergence  measured  in  the  energy  nor1,  as  d — >0.  More  precisely,  we 
def i ne 


n-model  we  will  understand  the  mlnimizer 


Bu(n) 


(B  (n) 
=  U1  ’ 


Bu<n), 


(2.7) 


e(d) 


A,  ,  B,B  (n), 
c  (u<j)  -  c  (  ud  ) 

cA(ud) 


1/2 


where  by  the  lower  index  d  we  express  the  dependence  of  the  solution  on  d. 
Assuming  B  =  A,  A  is  isotropic  Hooke’s  law  and  the  solution  Is  smooth,  then 


3 


does  not  lead  to  convergent  solutions  (a  =  0).  Nevertheless,  using  a 
modified  A,  i.e.  replacing  A  by  B.  we  get  the  value  a  =  1  again. 

This  modification  leads  to  the  well  known  Reissner-Mindl in  model  which  will 
be  discussed  in  Section  3. 

Let  us  consider  a  unit  square  plate  of  thickness  d. 

Cl  =  {  Xj.X2.X3  |  |Xj|  <0.5,  i  =  1,2,  Ix^l  <  ^  }  and  assume  that  the  plate 

is  simply  supported,  with  hard  support  and  is  uniformly  loaded.  We  assume 

that  the  plate  is  isotropic,  homogeneous,  with  a  Poisson’s  ratio  u  -  0.3 

7 

and  the  modulus  of  elasticity  E  =  10  .  The  solution  of  the  n-model  with  n 
small  is  smooth  and  essentially  has  no  boundary  layer  (in  contrast  to  the 
soft  simple  support). 


4 


In  Table  3.2  we  show  the  strain  energy  for  the  three  dimensional 
solution,  the  RM  model  (k  =  0.87)  and  the  (1,1,2)  model  (as  discussed 
later) . 

Table  3.2  The  convergence  of  the  RM  and  (1,1,2)  solutions 


d 

3DIM 

RM,  k  =-  0.87 

CRM 

(1,  1.2) 

C1  1  2 

0.  10 

0. 242115(-6) 

0.245521 (-6) 

11.8% 

0.240433 (-6) 

8.3% 

0.025 

0. 149115(-4) 

0. 149256 (-4) 

3.07% 

0. 149049(-4) 

2.  107. 

0.01 

0.232489 (-3) 

0. 232524( -3 ) 

I 

1 . 2% 

0.232472 (-3) 

0 . 85% 

By  e  and  e110  we  denoted  the  error  defined  by  (2.7)  for  the  RM  and  n 

KM  1  \c. 

(1,1,2)  model.  We  clearly  see  the  predicted  rate  of  convergence* 

For  more  detailed  analyses  of  various  boundary  conditions  we  refer  to  [2). 

3.  The  Relssner-Mlndl in  model  for  homogeneous  isotropic  materials. 

The  well-known  Relssner-Mindl in  model  is  described  by  the  following 
system  of  differentiable  equations  for  ^  and  u»  : 


D 

2 

( l-i>)A^1  + 

-  Kfid 

L  ♦  Js\ 

(3.1) 

D 

2 

( 1-v)A^2  + 

-  Kyd 

L  ♦ 

l*2  «*2J 

*cyd(Aw  +  ^-tf»)  +  q  =  C 

where 


D  = 


Ed' 


12(1  -  v ) 
E 


M  "  2(1  +  u) 


and  0  <  ic  S  1  is  a  shear  factor. 
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Here  0  denotes  the  rotation  vector  and  is  is  the  displacement  in  the 
vertical  direction. 

Further  we  define 


M 


11 


K22  = 


[3jh 


(3  2) 


*21  " 


=  M 


12 

°32 


D(s 
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[  dXi  3x2J 

n  &j>£\ 

2  u  [ax2  axj 

64?  *  *,) 
64?  *  *2} 


xpd 

xpd 


where  M^j,  i.j  =  1,2  are  the  moments  and  Q^,  J  =  1,2  are  the  shear 
forces. 

Let  n  and  s  be  the  unit  normal  and  right  oriented  tangential  vectors 

to  T  respectively.  Denote  further  0  =  n,  0  =  i^-s  and  define  M 

n  s  nn' 

Msg,  and  Mgn  according  to  (3.2)  replacing  x^  by  n  and  x^  by  s 
similarly  for  Q^n  and  Q^. 

At  the  boundary,  three  boundary  conditions  have  to  be  specified.  For 


example: 


i) 

Clamped  plate 

u %  - 

*  *z 

*  0 

ii) 

Soft  simple  support 

IS  = 

Mnn 

=  M  =  0 
ns 

ill) 

Hard  simple  support 

IS  = 

O* 

II 

M  =0 
nn 

iv) 

free 

M 

nn 

=  M 

ns 

0 

11 

c 

cP 

and  analogously  the  others. 

The  Reissner-Mindl in  (RM)  model  was  derived  using  various  mechanical 
principles.  It  has  been  shown  (see  e.g.  [3])  that  as  d — >0  the  solution 
of  the  RM  model  converges  to  the  solution  of  the  three-dimensional  model  in 
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the  (scaled)  energy  norm. 

Let  us  define  B  *  (b^>,  =  1..  .6  with 


V2"l 


*1+2^1  0 


V2“l 


where 


1  +  w  )(1-  2v: 


u  =  ,  Ei 

M1  “  2(1  +  v. 


V1  1  ♦  v 


„  _  E(1  ♦  2v) 

2 

1  (i  +i>r 


B  B 

and  denote  by  u  the  minimizer  of  G  (u)  (see  (2.3))  when  u  has  the  form 

(2.6)  with  ■  1,  n^  ®  0  and  H(C1)  is  defined  accordingly. 

For  example  in  the  case  of: 

i)  clamped  plate  :  K(fl)  =  (u  e  (hSq))2,  u  =  0  on  S> 

il)  soft  simple  support  :  X(Q)  =  (u  €  (H^UJ))2,  u^  =  0  on  S> 

iii)  hard  simple  support  :  H(£l)  =  (u  €  (H^Ul))2,  (u^.u^l’S  =  0, 

u^  =  0,  on  S> 

lv)  free  plate  :  Jf(fl)  =  (u  6  (H^(Q))2} 

B  B  B 

Let  now  <r  and  t  be  the  stresses  computed  from  u  by  using 
1  •  J  *  *  J 

Hooke’s  matrix  B  and  define 


d/2 

\i  '  'ii*: 

_  J  o 


'22*3d*3 
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(3.3) 


rc 

3M12  " 


T12X3dx3 


031  = 


032  =  T32dx3 


Now  It  is  noi.  hard  to  prove  the  following: 

Theorem  1.  We  have  with  (3.2)  and  (3.3) 


M11  "  M1J 


i  =  1,2 


031=  O3 j 


i  =  1,2 


B 

ut  -  u„ 


,  B 

*1  =  Ul, V 


where  M^,  Q}  ^  are  the  Relf  jner-Mlndl  in  moments  and  shear  forces  from 
(3.2).  The  proof  is  based  on  the  comparison  of  the  bil'near  forms  for  the 
variational  formulation  of  the  RM  model  and  the  (1,1,0)  model. 

We  remark  that  for  v  -  0,  the  matrices  A  and  B  are  identical. 


4.  The  Relssner-Mlndl in  model  for  homogeneous  orthotropic  plate. 


Let  the  Hooke’s  law  for  the  orthotropic  material  is 
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(4.1)  in  the  form  of  (2.1). 

Following  [4]  the  system  of  differential  equations  analogous  to 
equation  (3.1)  for  ^  and  m 


where 

C  =  {C^}  i.J.  =  1,2 

1  _v2\ 

Eu  E22 

u\2 

Eu  E22 


The  analog  to  (3.2)  is 
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(4.2) 
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and  M  ,  M  , 
nn  ns 
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ss 


^3n’  ^3s’  311,6  defined  accordingly 


B  = 


'11  12 


'21 

0 

0 

0 

0 


'22 

0 

0 

0 

0 


0 

0 

'33 

0 

0 

0 


0 

0 

0 

"12 

0 


0 

0 

0 

0 

kG. 


31 


0 

0 

0 

0 
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32  J 


Denoting  by  0^,  **  the  solution  of  the  Relssner-Mindl In  equation  for 

B  B 

given  boundary  conditions  and  by  u  the  minimlzer  of  G  ,  then,  as  in  the 
previous  section  the  Theorem  1  also  holds.  The  proof  is  analogous. 


5.  The  general  anisotropic  case. 

Consider  the  general  case  of  the  Hooke’s  law  defined  by  (2.1).  Then  if 
the  matrix  B  =  (b^}  is  such  that 

bu  -  *u  -  “S1  •  *-J  =  1-2-4 

b31  =  b32  =  b13  =  b23  =  b34  *  b43  =  0 
b33  *  a33 


biJ  =  KaiJ*  1 »  J  *  5,6 
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Then  once  more  Theorem  1  holds.  The  proof  Is  analogous. 


6.  The  h-p  version  of  the  finite  element  method. 

Let  the  domain  w  be  partitioned  into  two-dimensional  straight  or 


1 


in 


curvilinear  triangular  or  quadrilateral  (two-dimensional)  elements  e 

plane.  Then  the  three  dimensional  elements  are  e^  x  The»-p 

version  which  is  implemented  in  MSC/PROBE  and  STRIPE  program  uses  polynomials 
of  degree  q  (for  u^,  i  =  1,2,3)  in  the  x^  direction  and  polynomials  of 
degree  p  in  x^x2  plane  with  p  £  q.  For  quadrilateral  elements  both 
programs  use  the  serendipity  elements  in  the  xjx2  plane. 

The  elements  of  the  p.q  type  are  used.  Then  for  fixed  q  and  for 
p — »»  (or  when  the  size  of  elements  e^  converges  to  zero)  we  get  the 
solution  of  the  n-model  with  n  =  (q,q,q).  This  model  in  the  case  of 
symmetrical  load  is  equivalent  with  the  model  n  =  (n^.r^.n^),  n^  =  n2  = 


Tq+i' 

2 


1  and  n^  =  2 


where  by  [a]  we  denote  the  integral  part  of 
a.  It  Is  also  possible  to  use  different  q  In  different  elements.  Various 
material  properties,  isotropic,  anisotropic, are  available  in  the  programs. 
Hence  as  has  been  shown  in  previous  sections  the  RM  model  is  obtained  by 
properly  adjusting  the  Hooke’s  law  matrix. 

Let  us  consider  the  case  of  the  unit  square  with  hard  simple  support 
which  is  uniformly  loaded.  We  will  use  the  RM  and  (1,1,2)  model.  Let  us  be 
interested  li  the  locking  problem  l.e.  the  convergence  of  the  h-verslon 
using  x^x2  plane  elements  of  various  degrees. 

Consider  a  sequence  of  uniform  square  meshes  with  element  size  h.  In 
Figure  6. labc  we  show  the  error  (measured  in  the  energy  norm)  as  a  function  of 
h  for  various  h  and  various  degrees  of 

elements  p.  In  the  figure  we  also  depict  the  slope  of  the  asymptotic 
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relative  error  % 


error  rate  h 


Figure  6.1  The  convergence  of  the  h-version  for  the  RM  model,  a:  d  *  0.  1, 
b:  d  =  0.025,  c:  d  ■  0.01 

We  see  that  for  d  *  0. 1  the  rate  of  convergence  Is  a(p)  =  p.  For  d  * 
0.01  we  see  a(l)  *  0,  a(p)  <  p,  p  »  2,3,  but  a(4)  *  4  as  for  d  =  0.01. 
The  lower  rate  effect  Is  the  typical  locking  phenomenon.  For  d  ■  0.025  we 
see  locking  for  p  *  1,2  but  for  p  *  3  the  locking  is  essentially  not 
present . 
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Figure  6.2ab  shows  the  typical  locking  for  the  p-verslon.  We  see  here 
once  sore  that  for  Increasing  p  there  Is  no  difference  between  d  *  0.01 


and  d  *  0. 1  l.e.  that  the  p-verslon  is  locking  free. 


Figure  6.2  The  convergence  of  the  p-verslon  for  the  RM  model,  a:  d  *  0.1, 
b:  d  *  0.01 

Figure  6.3abc  shows  performance  of  the  h  version  for  the  model  (1,1,2) 
as  a  function  of  p,  analogously  as  before.  We  see  quite  analogous  results. 
Hence  for  p  >  4  the  method  Is  locking  free  also  here. 

In  [5]  we  have  addressed  the  general  model  of  locking.  We  have 
underlined  that  often  the  deterioration  of  the  finite  element  method  Is  the 
combined  effect  of  locking  and  decreasing  regularity  of  the  solution.  Hence 
it  Is  essential  to  guarantee  that  the  regularity  of  the  solutions  is 
essentially  uniform  with  respect  to  the  parameter  (in  our  case  d)  of 
interest.  We  have  for  this  reason  selected  the  hard  simple  support. 

In  the  case  of  the  model  (3,3,4)  the  solution  of  the  hard  simple 
support  Is  not  any  more  sufficiently  smooth  any  more  and  the  locking  Is  not 
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JtELATIVE  ERROR  % 


the  sain  factor  Influencing  the  convergence  rate.  This  leads  to  the  exactly 


1/4  1/8  1/12  1/16 


Figure  6.3  The  convergence  of  the  h-version  for  the  (1,1,2)  model,  a:  d  = 
0. 1.  b:  d  *  0.025  d  *  0.01. 


opposite  character  we  have  seen  above.  For  d  larger  we  see  smaller  rate  of 
convergence  than  for  d  smaller.  In  Figure  6.4  we  show  the  error  as  a 
function  of  h  for  p  ■  4.  It  clearly  Indicates  the  effect  of  the 
regularity  of  the  solutions  for  various  d.  The  effect  of  the  unsmoothness 
of  the  solution  for  d  *  0. 1  can  be  seen  by  comparing  the  accuracy  obtained 
for  different  meshes.  In  the  Table  6.  1  we  show  the  error  for  various  meshes 
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of  16  elements  for  d  ■  0. 1  and  d  ■  0.01.  These  aeshes  are  characterized  by 
coordinates  x^  i.  0. 

We  see  that  for  d  *  0. 1  the  uniform  mesh  (No.  1)  leads  to  a  larger 
error  while  a  refined  mesh  (No.  2)  leads  to  a  smaller  one.  On  the  other  hand 
the  uniform  aesh  (with  16  elements)  for  d  ■  0.01  is  the  optimal  one.  This 
clearly  explains  Figure  6.4  where  the  rate  of  convergence  for  d  *  0. 1  is 
smaller  (in  our  range)  than  for  d  *  0.01,  which  is  opposite  what  has  been 
seen  for  the  RM  and  (1,1,2)  model  before. 


—  h 


Figure  6.4  The  convergence  of  the  h-version  for  the  (3,3,4)  model. 


Table  6.1  The  energy  error  for  the  model  (3,3,4)  and  element  of  degree  4 
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COORDINATES 

EHR0R  % 

No 

■9 

xi 

X2 

X3 

n 

d  *  0.  1 

d  -  0.01 

B 

Bg| 

EO 

0.375 

0.50 

B 

0.209 

H 

ESI 

US 

0.45 

0.50 

HB 

0.434 

H 

Hi 

0.32 

0.43 

0.50 

— 

H 

0.30 

0.45 

0.49 

0.50 

B9 
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